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In eukaryotic cells, many motor proteins can move simultaneously on a single microtubule track. 
This leads to interesting collective phenomena like jamming. Recently we reported (P/ij/s. Rev. Lett. 
95, 118101 (2005)) a lattice-gas model which describes traffic of unconventional (single-headed) 
kinesins KIFIA. Here we generalize this model, introducing a novel interaction parameter c, to 
account for an interesting mechano-chemical process which has not been considered in any earlier 
model. We have been able to extract all the parameters of the model, except c, from experimentally 
measured quantities. In contrast to earlier models of intra-cellular molecular motor traffic, our model 
assigns distinct "chemical" (or, conformational) states to each kinesin to account for the hydrolysis 
of ATP, the chemical fuel of the motor. Our model makes experimentally testable theoretical 
predictions. We determine the phase diagram of the model in planes spanned by experimentally 
controllable parameters, namely, the concentrations of kinesins and ATP. Furthermore, the phase- 
separated regime is studied in some detail using analytical methods and simulations to determine 
e.g. the position of shocks. Comparison of our theoretical predictions with experimental results is 
expected to elucidate the nature of the mechano-chemical process captured by the parameter c. 



I. INTRODUCTION 

Motor proteins are responsible for intra-cellular trans- 
port of wide varieties of cargo from one location to an- 
other in eukaryotic cells [U, [1, 9 • One crucial feature of 
these motors is that these move on filamentary tracks (J| . 
Microtubules and filamentary actin are protein filaments 
which form part of a dual-purpose scaffolding called cy- 
toskeleton 5]; these filamentary proteins act like struts 
or girders for the cellular architecture and, at the same 
time, also serve as tracks for the intra-cellular transporta- 
tion networks. Kinesins and dyneins are two superfami- 
lies of motors that move on microtubule whereas myosins 
move on actin filaments. A common feature of all these 
molecular motors is that these perform mechanical work 
by converting some other form of input energy. How- 
ever, there are several crucial differences between these 
molecular motors and their macroscopic counterparts; 
the major differences arise from their negligibly small 
inertia. That's why the mechanisms of single molecu- 
lar motors [l], [3, H^and the details of the underlying 
mechano-chemistry [6[ have been investigated extensively 
over the last two decades. 

However, often a single filamentary track is used simul- 
taneously by many motors and, in such circumstances, 
the inter-motor interactions cannot be ignored. Funda- 
mental understanding of these collective physical phe- 
nomena may also expose the causes of motor-related dis- 
eases (e.g., Alzheimer's disease) thereby helping, pos- 
sibly, also in their control and cure. 

To our knowledge, the first attempt to understand ef- 



fects of steric interactions of motors was made in the con- 
text of ribosome traffic on a single mRNA strand 8] . This 
led to the model which is now generally referred to as the 
totally asymmetric simple exclusion process (TASEP); 
this is one of the simplest models of non-equilibrium sys- 
tems of interacting driven particles Ha [H. In the 
TASEP a particle can hop forward to the next lattice 
site, with a probability q per time step, if and only if the 
target site is empty; updating is done throughout either 
in parallel or in the random-sequential manner. 

Some of the most recent generic theoretical models of 
interacting cytoskeletal molecular motors [l^, [3, EBl 
are appropriate extensions of TASEP. In those models 
the motor is represented by a self-driven particle and the 
dynamics of the model is essentially an extension of that 
of the TASEP ^,Tr| that includes Lan gmuir-likc kinetics 
of attachment and detachment of the motors. Two differ- 
ent approaches have been suggested. In the approach fol- 
lowed by Parmeggiani, Franosch and Frey (PFF model) 
[13I . [l^ . attachment and detachment of the motors is 
modelled, effectively, as particle creation and annihila- 
tion, respectively, on the track; the diffusive motion of 
the motors in the surrounding fluid medium is not de- 
scribed explicitly. In contrast, in the alternative formu- 
lation suggested by Lipowsky and co-workers [H, |T3 , the 
diffusion of motors in the cell is also modelled explicitly. 

In reality, a motor protein is not a mere particle, but 
an enzyme whose mechanical movement is coupled with 
its biochemical cycle. In a recent letter [31 we consid- 
ered specifically the single-headed kinesin motor, KIFIA 
[H, do, [m, il m. The movement of a single KIFIA 
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motor had already been modelled earlier [20|, |2j| by a 
Brownian ratchet mechanism [25l [2^ . In contrast to 
the earlier models [13, [13, IT3. Il5| of molecular motor 
traffic, which take into account only the mutual inter- 
actions of the motors, our model explicitly incorporates 
also this Brownian ratchet mechanism of the individual 
KIFIA motors, including its biochemical cycle that in- 
volves adenosine triphosphate (ATP) hydrolysis. 

The TASEP-like models predict the occurrence of 
shocks. But since most of the bio-chemistry is captured 
in these models through a single effective hopping rate, it 
is difficult to make direct quantitative comparison with 
experimental data which depend on such chemical pro- 
cesses. In contrast, the model we proposed in ref. [3l 
incorporates the essential steps in the biochemical pro- 
cesses of KIFIA as well as their mutual interactions and 
involves parameters that have one-to-one correspondence 
with experimentally controllable quantities. 

Here, we present not only more details of our earlier 
calculations but also many new results on the properties 
of single KIFIA motors as well as their collective spatio- 
temporal organization. Moreover, here we also general- 
ize our model to account for a novel mechano-chemical 
process which has not received any attention so far in 
the literature. More specifically, two extreme limits of 
this generalized version of the model correspond to two 
different plausible scenarios of ADP release by the motor 
enzymes. To our knowledge j27| . at present, it is not pos- 
sible to rule out either of these two scenarios on the basis 
of the available empirical data. However, our new gen- 
eralized model helps in prescribing clear quantitative in- 
dicators of these two mutually exclusive scenarios; use of 
these indicators in future experiments may help in iden- 
tifying the true scenario. 

An important feature of the collective spatio-temporal 
organization of motors is the occurance of a shock or 
domain wall, which is essentially the interface between 
the low-density and high-density regions. We focus on 
the dependence of the position of the domain wall on 
the experimentally controllable parameters of the model. 
Moreover, we make comparisons between our model in 
the low-density regime with some earlier models of single 
motors. We also compare and contrast the basic features 
of the collective organization in our model with those 
observed in the earlier generic models of molecular motor 
traffic. 



II. DISCRETIZED BROWNIAN RATCHET 
MODEL FOR KIFIA: GENERAL FORMULATION 

Through a series of in-vitro experiments, Okada, Hi- 
rokawa and co-workers established [H, [2^, [2l|, [23, f23j 
that: 

i. KIFIA molecule is an enzyme (catalyst) and 
in each enzymatic cycle it hydrolyzes one ATP 
molecule; the products of hydrolysis being adeno- 
sine diphosphate (ADP) and inorganic phosphate. 
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FIG. 1: A biochemical cycle of a KIFIA motor (upper) and 
a three- valued discrete model for traffic of interacting KIFIA 
motors on a finite microtubule filament (lower). The states 
left to the dotted line in the upper figure correspond to 
strongly bound to microtubule states (state 1) while those 
right are weakly bound (state 2). denotes an empty site, 
and only 2 can move either to the forward or backward site. 
Transition from 1 to 2 occurs at the same site which corre- 
sponds hydrolysis, and the detachment also happens in this 
process. The attachment is possible only at the empty sites. 
At the minus and plus ends the probabilities are different from 
those at sites in the bulk. 



Thus, each biochemical cycle of a KIFIA motor 
consists of four states: bare kinesin (K), kinesin 
bound with ATP (KT), kinesin bound with ADP 
and phosphate (KDP) and, finally, kinesin bound 
with only ADP (KD) after releasing phosphate 
(Fig.©. 

ii. When a single-headed kinesin binds with a ATP 
molecule, its binding with its microtubule track is 
weakened by the ATP hydrolysis. Both K and KT 
bind strongly to microtubules. Hydrolysis of ATP 
leads to the state KDP which has a very short life- 
time and soon yields KD by releasing phosphate. 
KD binds weakly to a microtubule. After releasing 
all the products of hydrolysis (i.e., ADP and phos- 
phate), the motor again binds strongly with the 
nearest binding site on the microtubule and thereby 
returns to the state K. 

iii. In the state KD, the motor remains tethered to the 
microtubule filament by the electrostatic attraction 
between the positively charged if-loop of the mo- 
tor and the negatively charged ii^-hook of the mi- 
crotubule filament. Because of this tethering in the 
weakly bound state, a KIFIA cannot wander far 
away from the microtubule, but can execute (essen- 
tially one-dimensional) diffusive motion parallel to 
the microtubule filament. However, in the strongly 
bound state, the KIFIA motor cannot execute dif- 
fusive excursions away from the binding site on the 
microtubule. 

These experimental results for the biochemical cycle 
of KIFIA motors indicate that a simplified description 
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in terms of a 2-state model could be sufficient to un- 
derstand the collective transport properties. As shown 
in Fig. [1] one distinguishes a state where the motor is 
strongly bound to the microtubule (state 1) and a state 
where it is weakly bound (state 2). It is worth pointing 
out that such a simplified 2-state model, however, may 
not be adequate to capture the biochemical cycle of other 
motors like, for example, conventional kinesins. In such 
situations, a more detailed 4-state model is required. 

As in the TASEP-type approach of the PFF model, 
the periodic array of the binding sites for KIFlA on 
the microtubule are represented as a one-dimensional 
lattice of sites that are labelled by the integer index i 
[i — 1, L). KIFlA motors are represented by particles 
that can be in two different states 1 and 2, corresponding 
to the strongly-bound and weakly bound states. To 
account for the empirical observations, the model also 
contains elements of a Brownian ratchet. As in the 
PFF model, attachment and detachment of a motor are 
modelled as, effectively, creation and annihilation of the 
particles on the lattice. We use the random sequential 
update, and the dynamics of the system is given by the 
following rules of time evolution: 

(1) Bulk dynamics 

If the chosen site on the microtubule is empty, i.e., in 
state 0, then with probability ojadt a motor binds with 
the site causing a transition of the state of the binding 
site from to 1. However, if the binding site is in state 
1, then it becomes 2 with the probability ojhdt due to 
hydrolysis, or becomes with probability ujddt due to 
the detachment from the microtubule during hydrolysis. 

If the chosen site is in state 2, then the motor bound 
to this site steps forward to the next binding site in front 
by a ratchet mechanism with the rate ojfov stays at the 
current location with the rate ujs- Both processes are 
triggered by the release of ADP. How should one modify 
these update rules if the next binding site in front is 
already occupied by another motor? Does the release of 
ADP from the motor, and its subsequent re-binding with 
the filamentary track, depend on the state of occupation 
of the next binding site in front of it? To our knowledge, 
experimental data available at present in the literature 
are inadequate to answer this question. Nevertheless, we 
can think of the two following plausible scenarios: in the 
cases • • • 21 • • • , or • • • 22 • • • , the following kinesin, which 
is in state 2, can return to state 1, only at its current 
location, with rate lOs if ADP release is regulated by the 
motor at the next site in front of it. But, if ADP release 
by the kinesin is independent of the occupation status of 
the front site, then state 2 can return to state 1 at the 
fixed rate LOs+ujf^ irrespective of whether or not the front 
site is occupied. 

Therefore, we propose a generalization of our origi- 
nal model by incorporating both these possible scenarios 
within a single model by introducing an interpolating pa- 
rameter c with < c < 1. In this generalized version of 
our model, a motor in the state 2 returns to the state 1 



at the rate + (1 — c)a;/. The parameter c (0 < c < 1) 
allows interpolation between the two above mentioned 
scenarios of ADP release by the kinesin. For c = 1 the 
transition from the strongly to the weakly bound state in 
the ratchet mechanism depends on the occupation of the 
front site. This is the case that has been treated in p^ . 
where the release of ADP by a nucleotide-bound kinesin 
is tightly controlled by the kinesin at the next binding 
site in front of it. On the other hand, for c < 1 the tran- 
sition rate will depend partially on the occupation of the 
front site. For c = the ADP release process becomes 
completely independent of the state of the preceeding 
site. 

As long as the motor does not release ADP, it executes 
random Brownian motion with the rate cjf,. 

(2) Dynamics at the ends 

The probabilities of detachment and attachment at the 
two ends of the microtubule can be different from those 
at any other site in the bulk. We choose a and 5, instead 
of uJa^ as the probabilities of attachment at the left and 
right ends. Similarly, we take 71 and /3i, instead of ojd, 
as probabilities of detachments at the left and right 
ends, respectively (Fig. [T]). Finally, 72 and /32, instead 
of LJb, are the probabilities of exit of the motors through 
the two ends by random Brownian movements. 

For the dynamical evolution of the system, one of the 
L sites is picked up randomly and updated according to 
the rules given below together with the corresponding 
probabilities (Fig. [2|): 



Attachment : 
Detachment 
Hydrolysis : 



Brownian motion 



Ratchet 




■ 1 with LOadt 

with ojddt 

2 with ojhdt 

20 ^ 02 with bJbdt 
02 20 with LObdt 

10 with LOadt 

IX with {uja + {l-c)ujf)dt{b) 

01 with LOfdt 



(1) 
(2) 
(3) 

(4) 



Here X denotes an occupied site irrespective of the chem- 
ical state of the motor, i.e., a site occupied by a motor 
that is in either state 1 or state 2. 

The ratchet mechanism (O is triggered by the release 
of ADP and summarizes the transitions of a particle from 
state 2 to state 1. It distinguishes the two initial states 
20, where the front site is empty, and 2X^ where the 
front site is occupied. We see that the overall transition 
rate from state 2 to state 1 is -t- w/ if the front site 
is empty (initial state 20), and it is cj^ + (1 — c)LOf if 
the front site is occupied (initial state 2X). This reflects 
the dependence of the ADP release rate on the front site 
occupation whenever c ^ 0. 

The physical processes captured by the rate constants 
ujf and LOs can be understood as follows by analyzing 
the Brownian ratchet mechanism illustrated in Fig. O 
For the sake of simplicity, we consider only one molecu- 
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FIG. 2: Schematic description of the three-state model of a 
single-headed kinesin motor that follows a Brownian ratchet 
mechanism. In the special case 2X — > IX, which has not been 
shown explicitly for the sake of simplicity, the rate constants 
would get modified following the prescriptions described in 
the text. 
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FIG. 3: The two forms of the time-dependent potential used 
for implementing the Brownian ratchet mechanism. 



lar motor, and let us imagine that the potential seen by 
the motor periodically oscillates between the sawtooth 
shape and the flat shape shown in Fig. [31 When the 
sawtooth form remains "on" for some time, the particle 
settles at the bottom of a well. Then, when the potential 
is switched "off" , the probability distribution of the po- 
sition of the particle is given by a delta function which, 
because of free diffusion in the absence of any force, be- 
gins to spread. After some time the Gaussian profile 
spreads to such an extent that it has some overlap also 
with the well in front, in addition to the overlap it has 
with the original well. At that stage, when the sawtooth 
potential is again switched on, there is a non-vanishing 
probability that the particle will find itself in the well 
in front; this probability is proportional to the area of 
the hatched part of the Gaussian profile shown in Fig. [3] 
and is accounted for in our model by the parameter tof. 
There is also significant probability that the particle will 
fall back into the original well; this is captured in our 
model by the parameter cOs- 



The detachment rate Ud — 0.1 s~^ is found to be inde- 
pendent of the kinesin population. On the other hand, 
uja = 10^ C/M-s depends on the concentration C (in 
M) of the kinesin motors. In typical eucaryotic cells 
in-vivo the kinesin concentration can vary between 10 
and 1000 nM. Therefore, the allowed range of uja is 0.1 
s"^ < Wa < 10 s-\ 

Total time taken for the hydrolysis of one ATP 
molecule is about 9 ms of which 4 ms is spent in the 
state 1 and 5 ms in the state 2. The corresponding 
rates 1/4 and 1/5 are shown in Fig. [H The motion of 
KIFlA is purely diffusive only when it is in the state 2 
and the corresponding diffusion coefficient is denoted by 
the symbol D2. Using the measured diffusion constant 
D = 40,000 nmVs [20| and the relation D2 = (9/5)1), 
we get 132 = 72, 000 nmVs (see Fig.gjb)). The time 
must be such that Wf, ~ Z?2/(8nm)^, and, hence, we get 
ujb ^ 1125 s-^ 

Moreover, from the experimental observations that the 
mean step size is 3 nm whereas the separation between 
the successive binding sites on a microtubule is 8 nm, we 
conclude luj /us — 3/8. Furthermore, from the measured 
total time of each cycle, we estimate that lUs + ujf — 
200 s~^. From these two relations between w/ and usg 
we get the individual estimates cOs — 145 and i^f — 
55 s-\ 

Assuming the validity of the Michaelis-Menten type 
kinetics for the hydrolysis of ATP [5] , the experimental 
data suggest that 
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where [^TP] is the ATP concentration (in mM), Km is 
the Michaelis constant given by Km = 0.1 mM in this 
case. V and Vmax (in ms~^) are the reaction rate and 
its maximum value respectively. As mentioned earlier 
^/Vmax — 9 ms. Since 1/V — uj^^ + 5 ms, we finally get 



0.1 mM 



ATP concentration (in mM) 



ms (7) 



so that the allowed biologically relevant range of ujh is 
0<LUh< 250 s~i. 

Up to now, experimental investigations could not de- 
termine the parameter c. We therefore treat it as a free 
parameter in the following to study the effects that it 
has on the phase diagram, position of shocks etc. Com- 
parison with empirical results then might help to get an 
estimate for c. 



B. Mean- field equations 



Parameters of the model 



From experimental data [T^, [2l|, good estimates for 
the parameters of the suggested model can be obtained. 



Let us denote the probabilities of finding a KIFlA 
molecule in the states 1 and 2 at the lattice site i at time 
t by the symbols Si and Wi, respectively. In mean-field 
approximation, the master equations for the dynamics of 
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FIG. 4: The biochemical cycle of KIFIA is shown to define 
some important parameters which can be extracted from ex- 
perimental data. See text for more details. 



the interacting KIFIA motors in the bulk of the system 
are given by 

^ = LOa{'^- S.i-Wi)~UJhSi-UJdSi 

at 



dt 



+ {l-c)ujfW,{S,+i + W,+i), 



(8) 



-{ws + (1 - c)ujf)W,{S,+i + W,+i) 
-i^bWd^ - S,+i - W,+i ~ - iy,_i) 
+u;b{W,.i+W,+i)il-S,-W,) 

-UbW,i2 - S,+i - W,+i - ^,„i - W,^i) 
+iObiW,-i+W,+i)il~S,-W,). (9) 

The corresponding equations for the left boundary (i = 
1) are given by 



- ail-Si~Wi)+LJ,Wi~LJhSi-^iSi 



+ {l-c)ujfWi{S2+W2), 



(10) 



dWi 
dt 



LUhSi - {lUs + UJf)Wi + CUfWi{S2 + W2) 
--f2Wi+LObW2[l- Si-Wi) 

-ujbWiil- S2-W2), (11) 



while those for the right boundary (i = L) are given by 
^ = 5{\- SL-WL)+iOfWL-i{l- Sl-Wl) 



dWL 
dt 



+UJsWl - UJhSl - PiSl, 
= LOhSL - UJsWl - P2WL 

+iUbWL-lil~ Sl~Wl) 
-UJbWLil-SL-l-WL^l). 



(12) 



(13) 



In the following we shall determine solutions of this set 
of equations for several cases and compare with the cor- 
responding numerical results from computer simulations. 



III. COMPARISON WITH OTHER MODELS 
FOR MOTOR TRAFFIC 

In this section we compare our model with earlier mod- 
els of molecular motor traffic. The first two subsections 
describe models developed for non-interacting molecu- 
lar motors whereas in the last subsection we collect the 
main results for the PFF model which has been intro- 
duced to study collective effects in motor traffic. A more 
detailed comparison with models of interacting motors 
will be taken up later in section |V] of this paper. 

Chen [2^ developed a model for single-headed kinesins 
assuming a power stroke mechanism. He assumed that 
each kinesin can attain three distinct states which were 
labelled by the symbols 0, 1 and 2. The kinesin was 
assumed to be detached from the microtubule in the state 
0, but boimd to microtubule in the other two states. The 
states 1 and 2 were assumed to differ from each other 
by the amount of their tilt in the direction of motion. 
The molecule steps ahead by exactly 8 nm in one cycle 
consuming one ATP molecule. This power-stroke model 
fails to account for several aspects of experimental data 
(for example, the distribution of the steps sizes, including 
backward steps) on KIFIA and, therefore, will not be 
considered further for quantitative comparison. 



Comparison with Sasaki's Brownian ratchet 
model 



In contrast to t he p ower-stroke model developed by 
Chen [2^, Sasaki [24| quantified the Brownian-ratchet 
model for a single KIFIA motor proposed by Okada and 
Hirokawa Il9l [201 . He used the standard Fokker-Planck 
approach [25ll26]. In this formulation, the particle, which 
represents a kinesin, is assumed to be subjected to a time- 
dependent periodic potential as given in Fig. [31 The po- 
tential switches from one shape Vi{x) to another shape 
V2{x) with rate wi and the reverse switching takes place 
at a rate UJ2 ■ One of the shapes of this potential Vi (x) 
is taken to be a periodic repetition of a saw-tooth where 
each saw-tooth itself is asymmetric. Suppose, the height 
of the maximum of each sawtooth is U. The shape of the 
form of the potential V2(x) was assumed to be flat, i.e., 
V2{x) — for all X. Sasaki calculated the average speed v 
and the diffusion coefficient D as functions of U, uji and 

1^2- 

One advantage of our model over Sasaki's model is 
that we do not make any ad-hoc assumption regarding 
the shape of the potential as the potential does not enter 
explicitly into our formulation. It is possible to identify 
LUi in Sasaki's model with coh in our model. The rate 
constant UJ2 can be related to the rates in our model in 
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the following way: uj2 = + oj/ if the preceding site is 
unoccupied and = + (1 ^ '^)^f occupied. 



B. Comparison with Fisher-Kolomeisky multi-step 
chemical kinetic model 



Next we make a comparison between our model and 
the multi-step chemical kinetic approach developed by 
Fisher and Kolomeisky (29l. fsol. Isij] for molecular motors. 
In the simplest case of a single filament, the equispaced 
binding sites on a microtubule are assumed to form a 
one-dimensional lattice. It is assumed that there are M 
distinct discrete intermediate chemical states on a bio- 
chemical pathway between two consecutive binding sites. 
The motor in state ji (i.e., in chemical state j located at 
spatial position i where l<j<M,l<i<L) can make 
transitions to the states + and with the rates 

Uj and Wj^ respectively (see Fig. [5])- Note that we have 
labelled the chemical states in such a way that Mi = li+i 
[Mi = 2 in Fig.O such that, completion of the chain in 
forward (backward) transitions through these M states 
would translocate the motor forward (backward) by one 
lattice spacing. 



i-1 




FIG. 5: Fisher-Kolomeisky multi-step chemical kinetic model 
of molecular motors. 

Clearly, in the absence of attachment and detachment 
of the motors, our model for a single KIFlA reduces to 
the Fisher-Kolomeisky multi-step chemical kinetic model 
of molecular motors on a single filament (see Figl6]) where 
M = 2, as emphasized by a slight redrawing of our model 
in Fig. H 

Direct quantitative comparison with our model is also 
possible. For example, in the special case where only 
forward transitions are allowed and M — 2, the average 
speed of the motor in the Fisher-Kolomeisky model is 
given by 



UiU2 
Ul + U2 



(14) 



where distance is measured in the units of spacing be- 
tween two successive binding sites (8 nm in case of mi- 
crotubule). In Sec. lTVBl we will derive an analogous ex- 
pression for our model, see ([TO]) . 



i+I 






FIG. 6: In the absence of attachment and detachment, our 
model is equivalent to Fisher-Kolomeisky model shown in 
Fig. El 



C. Comparison with PFF-model 

The Parmeggiani-Franosch-Frey model (PFF model) 
combines the TASEP with Langmuir kinetics. The 
motors are assumed to step forward one site with rate 
p if the front site is empty, but do not move if this site 
is occupied (exclusion). A backwards movement is not 
possible. In addition, motors can attach to empty sites 
with rate uja and detach from a site with rate LOd- This 
might be the simpliest model for intracellular transport 
including adsorption and desorption. Although quite ba- 
sic, it already reproduces the qualitative behavior of a 
large class of many-motor systems. It not only shows 
high-, low- and maximum-current phases like TASEP, 
but also phase coexistence for distinct parameter ranges, 
while phase domains are separated by stationary domain 
walls (shocks). These shocks are also observed in ex- 
periments |18j. Shock phases appear if the Langmuir 
kinetics are of the same order as motor attachment and 
detachment at the ends. It means that in the continious 
limit where system size i — > oo, the local attachment 
and detachment rates coa and ujd have to be rescaled so 
that the global attachment and detachment rates defined 
as fta ■= i^aL, fid ■= oJdL stay constant. One can argue 
that the topology of the phase diagram of the PFF model 
is quite universal for systems that, as the PFF model, 
possess a current-density relation with one single maxi- 
mum and the same Langmuir kinetics [isj, so even more 
complex models might show similiar qualitative behavior 
as the PFF model. 

Although the PFF model reproduces qualitative prop- 
erties of intracellular transport quite well, it is difficult 
to associate the hopping parameter p quantitatively with 
experimentally accessible biochemical quantities because 
the biochemical processes of a motor making one step 
are usually quite complex. The PFF model does not 
take into account these processes. Furthermore, it is not 
possible to include interactions in the PFF model that 
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ATP (mM) 


L^h (1/s) 


V (nm/ms) 


D/v (nm) 


r (s) 


oo 


250 


0.201 


184.8 


7.22 


0.9 


200 


0.176 


179.1 


6.94 


0.3375 


150 


0.153 


188.2 


6.98 


0.15 


100 


0.124 


178.7 


6.62 



TABLE I: Predicted transport properties in the low-density 
limit for four different ATP densities, r is calculated by av- 
eraging the intervals between attachment and detachment of 
each KIFIA. 



only influence particular transitions of the biochemical 
states of the motor. The advantage of our model is the 
possibility of calibration of the model parameters with 
experimentally controllable parameters ATP- or motor 
protein concentration. Through the parameter c we can 
include, at least phenomenologically, an interaction that 
controls the transition from one state (2) to another (1). 



IV. SINGLE MOTOR PROPERTIES AND 
CALIBRATION 



In this section we first investigate the dynamics of our 
model in the limit of vanishing inter-motor interactions. 
This helps us to calibrate the model properly by com- 
paring with empirical results. Then we compare the 
non-interacting limit of our model as well as the cor- 
responding results with earlier models of non-interacting 
motors to elucidate the similarities and differences be- 
tween them. 



Calibration of our model in the low-density 
limit 



B. Non-interacting limit of our model: a 
mean-field analysis 

For the case of a single KIFIA molecule, all interaction 
terms can be neglected and the mean- field equations (O, 
([9]) for the bulk dynamics are linearized and simplify to 



dSj 

<m 

dt 



= OJa{l -Si- W.i) + OJfWi-l + WsWi 

- L^hSi - uJdSi, (15) 

= LUhSi - UJsWt - LOfWi 

+ L0b{Wi.i+W,+i)-2L0bW,. (16) 



The boundary equations P^ - (fTT|) also get simplified in 
a similar way. 

Assuming periodic boundary conditions, the (homo- 
geneous) solutions {Si.Wi) = {S,W) of the mean-field 
equations ([TS]), in the steady-state are found to be 



S = 
w = 



UJai^^s + W/) 



t^gUJh 

LOaiuJh + t^s + ^/) + ^d(t^s + ^/) ' 

The corresponding flux is given by 

J = LUfW 

LUf{uJa + ^d) + ^a{^s + ^/i) + ^d^s 
_ iOh 

[1 + K) + {n, + nh) + n,K' 



(17) 
(18) 



(19) 



where K = ujd/oJa, = ^h/^f and fig = lOsl^j- 

If we make the correspondence u\ — ujh, and U2 = (^f 
the expression (|19p for the average speed of KIFIA in 
our model, in the special case lUs — (i.e., no reverse 
transition) reduces to the Fisher-Kolomeisky result (|14p . 
A more general version of the non-interacting limit of our 
model is treated in appendix A. 



An important test of our model would be to check if 
it reproduces the single molecule properties in the limit 
of extremely low density of the motors. We have already 
explained earlier how we extracted the numerical values 
of the various parameters involved in our model. The 
parameter values Ua = ct — 1.0 x 10~^ s~^, allows real- 
ization of the condition of low density of kinesins. Using 
those parameters sets, we carried out computer simula- 
tions with microtubules of fixed length L — 600 which is 
the typical number of binding sites along a microtubule 
filament. Each run of our simulation corresponds to a 
duration of 1 minute of real time if each timestep is inter- 
preted to to correspond to 1 ms. The numerical results of 
our simulations of the model in this limit, including their 
trend of variation with the model parameters, are in ex- 
cellent agreement with the corresponding experimental 
results (see Table |T|. 



V. COLLECTIVE FLOW PROPERTIES 

In the following we will study the effects of interac- 
tions between motors which lead to interesting collective 
phenomena. 

A. Collective properties for c — 1 

We first look at the case c = 1 originally studied in 
[TH - In mean- field approximation the master equations 
([5]), ©for the dynamics of the interacting KIFIA motors 
in the bulk of the system are nonlinear. Note that each 
term containing w/ is now multiplied by the factor of 
the form {1 — Si — Wi) which incorporates the effects of 
mutual exclusion. 
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Assuming periodic boundary conditions, the solutions 
{Si,Wi) = iS,W) of the mean-field equations dH), ^ in 
the steady-state for c = 1 are found to be 



S 
W 



^n,, - - {fls - 1)K + ^fD 



2K{1 + K) 

+ {Qs + ^)K-^/D 



2K 



where K = ujd/^a, = ^^h/^f, = <^sl^f-, and 



(20) 
(21) 

(22) 



Thus, the density of the motors, irrespective of the in- 
ternal "chemical" state, attached to the microtubule is 
given by 



p = S + W 



+ ns + i^s + i)K - Vd + 2 

2{1 + K) 



(23) 



This is the analogue of the Langmuir density for this 
model; it is determined by the three parameters K, il/j 
and fig. Note that, as expected on physical grounds, 
S + W ^ 1 as K ^ whereas S + W^OasK^oo. 
The probability of finding an empty binding site on a 
microtubule is KS as the stationary solution satisfies the 
equation S + W + KS ^ 1. 

The steady-state flux of the motors along their micro- 
tubule tracks is given by 

J ^ujfW{l- S -W). (24) 

Using the expressions (PT|) for S and W in equation 
for the flux we get the analytical expression 



J = 



4Kil + K) 



(25) 



FIG. 7: Fundamental diagram (i.e., flux-versus-density rela- 
tion) for the traffic flow of KIFIA in our model. 



The flux obtained from the expression (|25p for several 
different values of tUh are plotted as the fundamental dia- 
grams for this model in Fig. [71 Note that, in general, this 
model lacks the particle-hole symmetry. This is obvious 
from the flux can be recast in general as 



J = 



UJh+UJs+ t^/(l - Cp) 



(26) 



This is easily derived by substituting the relation p 
S + W and the constant solution of © 



{ujs+i^f)W + LUh{p-W) + cujfWp^O (27) 



into the definition of the flux 



Next we consider two limiting cases. In case I (w/ ^ 
LUh — ujs) the forward movement is the rate-limiting pro- 
cess and in case II {cuh cof — ujs) the availability of 
ATP and/or rate of hydrolysis is the rate-limiting pro- 




1. Case I (ujf <^ LUh — ^s) 



In this case, 



S 



1 



2K{1 + K) 



1 



1 



if 4 



risil + K) + K+-il + K){3 + K) (n, - K') - -17,(1 + KY{i + K)' + K' - — 
Z a 2\lh 



(28) 
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W 



1 
2K 



ns{i + K) 



K--il 



+ K)i3 + K) [VLs - K^) + -0,(1 + Kf{2, + Kf - 



(29) 



so that the total density is 

P = 5 + M..^±^M1±4I±4_(V^^. (30) 

2{1 + K) ^ ' 

Therefore, in this case, the steady-state flux is given by 

^^.,T^(l-p)^-^[^(^ + ^^-^](^'^). (31) 

K 

In this case, in addition, if ii' = ^ ^ 1, i.e., detachments 
are rare compared to attachments, can be treated as 
negligibly small and, hence, equations (p8)) . (p9)) and pO)) 
simplify to the forms 



1 



2(1 + K) 



(32) 



2. Case 11 (uju cof lOs) 



In this case also the flux (|26p can be interpreted to 
be that for a TASEP where the particles hop with the 
effective effective hopping probability 



(2) 

9cff 



(39) 



that depends on the density p. The specific form of (/off 
in equation (|39|) is easy to interpret physically. A tightly 
bound motor attains the state 2 with the rate ujh and 
only a fraction — , — r of all the transitions from the 
state 2 lead to forward hopping of the motor. 



B. Collective properties for c — 



w 



K 



The corresponding formula for the flux becomes 
J ^ - P) 



where 



(1) ^ + _ ^ 
2 + K 2 



(33) 



(34) 



(35) 



(36) 



Note that this effective hopping probability is also de- 
rived directly from ^2E\i by putting ujf -^uju — . 

Thus, the result for the flux in the special case can 
be interpreted to be that of a system of "particles" hop- 
ping from one binding site to the next with the effective 
hopping probability q^^^ . 

However, if we assume only ujf, but the relative 

magnitudes of LOh and lo^ remains arbitrary, 



fis - {iis - i)K + nh + n^ii + K)-K 



2K(l + K) 



nh + ns + + i)K - Qh - (1 + K)n, + k 

2K 



(37) 

""11 I ""s I I ^y^' ""11 ' -"y-s I ' ' 
W ~ — 1 • 

(38) 

Physically, this situation arises from the fact that, be- 
cause of fast hydrolysis, the motors make practically in- 
stantaneous transition to the weakly bound state but, 
then, remain stuck in that state for a long time because 
of the extremely small rate of forward hopping. 



We now consider the case c — where ADP release 
by the kinesin is independent of the occupation status of 
the front site. Let us study the stationary state of the 
mean- field equations ([8]), ([9]) in the case c = 0. From ([9]) 
we get 



UJs 



UJf 



(40) 



by neglecting the terms that represent Brownian motion. 
Substituting this into ([8]) we have 

ujfH,.i{l-H,)-LOjH,{l-H,+i)-adH,+aa{l-Hi) = 0, 

(41) 

where we put 



UJs+LOf + UJh 

UJs + UJf 
UJh 

UJs + UJf + UJh 

L 

UJh 



UJd 



ad , 



(42) 
(43) 
(44) 



Eq. (|¥T|) is the same equation as for the stationary PFF 
model. Therefore, the phase diagram of this model would 
be identical to that of the PFF model in mean field ap- 
proximation if we rescale all the parameters by 
One has to stress that this model is not exactly identical 
to the PFF model. While mean field approximation is ex- 
act for the PFF model in the continuous limit, our model 
shows correlations [s^ that lead to different density pro- 
files and phase diagrams (see Sec. IVI A| . Nevertheless, 
the topological structure of the phase diagrams remains 
the same in both models and the differences are not quite 
large. 

So far we have discussed two possible scenarios of ADP 
release by kinesin; in one of these the process depends on 
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(ATP[mM]) 

0,2 .. 

(0.9) 




(KIFlA[nM]) 



0.00001 (1) 0.00005(5) 0.001 (100) 



FIG. 8: (Color online) Space-time plot of the model system 
for c — 1. Each row of squares represents the state of the 
system at one single instant of time whereas successive rows 
(in the upward direction) correspond to the state of the sys- 
tem with increasing time. The blue and red squares indicate 
kinesins in the states 1 and 2, respectively, while the white 
squares correspond to empty binding sites on the microtubule. 
Total number of binding sites is 600, and the configurations 
of the system are displayed for the last 1200 time steps of 
a simulation run up to a total of 2 x 10^ time steps, start- 
ing from an initial state where all the binding sites on the 
microtubule were empty. The other model parameters are 
LUa = 0.3,uJd = 0.2, ujh = 400,a;/ = 600,ujs = 200,tJb = 50 
for bulk, and a = 50,/3i = /32 = 700,7i = 72 = 5 = for 
boundaries. 



the status of occupation of the target site (c = 1) whereas 
it is autonomous in the other {c — 0). To our knowledge, 
at present, the available experimental data can not rule 
out either of these two scenarios of ATP hydrolysis by 
kinesins. Therefore, we have introduced the parameter 
c that interpolates both these possible scenarios. As we 
have seen in this section, the extended model interpo- 
lates, at least on the level of mean-field theory, between 
the PFF model and the model introduced in 18]. In the 
following section we will discuss some properties of the 
extended model including case < c < 1 in more de- 
tail. We focus on the density profiles and especially the 
properties of shocks. 



VI. POSITION OF THE SHOCK 

One of the interesting results of the model is the exis- 
tence of a domain wall separating the high-density and 
low-density phases in the steady state of the system. One 
such configuration is shown in the space-time diagram in 
Fig. m In this section we shall determine the position 
of the shock, i.e., the domain wall, and the trends of its 
variation with the model parameters Wa and u f , etc. 



A. Analytical treatment in the continuum limit 

Let us first introduce the variable x — -fjE^; since 1 < 
i < L, we have < a; < 1. We map our system {Si, Wi) 
into {S{x),W{x)), and consider the continuum limit by 
considering L to be large enough: 

' ^ * (45) 



S{x±e)^Six)±e— 
ox 



2 dx^ 

for S{x ± e) and a similar expansion for W{x ± e), where 
e = 1/L. Using this Taylor expansion, we get 

dS{x,t) 



dt 



Wa(l - 

u;f[W 



dx 

, ^dSix,t) _^^dW{x,t)^ 



dx 



dW{x,t) 



dt 



S + W 



avfW 

- {ujs+ujf)W + ujhS . 
In the stationary state, we have 
dS{x) Us+ujf LOh S 



dx 

dS{x,t) 



dx 



, dWix,t) 
e H e 



dx 



(46) 



dx 
dW{x) 
dx 



CUJf 

1 



CLUf W 



dx 



.(nfi- 

1 -S-Wuf 

UJh + c{uJa + UJd) 
CUJf 



UJf — CUJn 



CUJf 

S-WiS + W). 



-W 



(47) 



Moreover, from the left boundary equations, by letting 
Si = S2 = S{0) and Wi ^ W2 = W{0), we obtain 

a[l - S{0) - W{0)] - uJhS{0) + ujsW{0) 

+ {l~c)ujfW{0)[S{0) + W{0)]^0 (48) 

-{LJs+UJf)W{0)+UJhS{0) 

+cujfWiO){S{0) + W{0))^0, (49) 



and, hence. 



^(0) 



[ca{a - uis)/ujf] 
ca + uJh 



W{Q) = — . 

Similarly from the right boundary conditions 

-ujhS{l)+ujsW{l)-fiS{l) 

+ujfW{l)[l- S{l)-W{l)]=Q, 
-uj,W{l)+ujhS{l)-pW{l)^Q- 
Solving these equations we have 



(50) 



(51) 
(52) 



5(1) 



UJh 



u}h+ujs+ (3 



£_ 

UJh+UJs + (3 UJf 

UJf ■ 



(53) 
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I^a^O.Ol 


UJa=0.025 


UJa = 0.05 


LJa=0.065 


200 


0.725 


0.5 


0.318 


0.253 


150 


0.776 


0.571 


0.382 


0.311 


125 


0.808 


0.618 


0.425 


0.35 



TABLE II: The position of the locahzed shock. Parameters 
are L = 600, c = 1, uid = P = 0.1, a = uia, ijJf = 145 and 
ujs = 55. 

Note that the pair of coupled equations ([47|) involves 
only the first order derivatives of S and W with respect 
to X whereas we have two sets of boundary conditions 
and (|53p. Therefore, if we integrate the equations 
(I47p using the boundary conditions (jSOp , the solution may 
not, in general, match smoothly with the other solution 
obtained for the same equation using the boundary con- 
ditions ((53|) . The discontinuity corresponds to a shock or 
domain wall. 

The continuity condition gives J/ — Jr where the flow 
just at the left side is denoted by J/ = lu jhi{\~ri—hi) and 
that at the right is Jr- Thus we integrate (|47l) numerically 
by using ([50]) to the right end, and we also integrate them 
by (|53p to the left end, and seek the point where J; = Jr 
is attained. 

We have investigated the shock position by changing 
the values of u>a and Uh- The results are given in Ta- 
ble [ni in the case c = 1 , which quantitatively agree with 
numerical simulations as shown in Fig. (S] 



Simulation 

Mean Field theory 



0.8 - 



0.6 - 



0.4 - 



0.2 - 




FIG. 9: (Color online) An example of a density profile with 
shock obtained by integrating the mean-field equations as well 
as that of numerical result. Parameters are uih ~ 200, tOa ~ 
0.01. 



B. Shock position from simulations 

In this subsection we locate the position of the shock 
in our model using a new shock tracking probe (STP) 
which is an extension of "second class particles" (SCP) 



[3J| used earlier for locating domain walls in computer 
simulations of driven-diffusive lattice gas models defined 
on a discrete lattice. In the standard TASEP model, a 
SCP is defined as one that behaves as a particle while 
exchanging position with a hole and behaves as a hole 
while exchanging position with a particle. As a result, 
the second class particle has a tendency to get localized 
at the domain wall (or, the shock). Other types of STP 
have also been considered in the literature [35] 

The rules for the movements of the STP in our model 
of KIFlA traffic have been prescribed by extending those 
for SCP in TASEP. Let us use the symbols 1 and 2 to 
denote the STPs which correspond to the states 1 and 
2, respectively, of the particles. Now, in the special case 
c = 1 , we define the following rules for the movements of 
the STPs: 



1 


2 


with rate LUh 


2 ^ 


1 


with rate lUs 


20 - 


01 


with rate w/ 


20 - 


02 


with rate Wf, 


02 - 


20 


with rate tUb 


OX 


■■XX ~ 


~> XX ■ ■ ■ XX with rate oja 


2X 


XI 


with rate ujf 


2X 


-> X2 


with rate u>b 


X2 


2X 


with rate ojb 



with X and X denoting occupation in either state of 
particles or STPs respectively, while • • • denotes a line 
of sites occupied by particles. Further extension of these 
rules for arbitrary c is straightforward. 

These rules satisfy the STP-principle: if the selected 
site is a STP it behaves like a particle, while if the se- 
lected site is a particle it treats STPs in its vicinity as 
holes (by changing sites respectively). Note that there 
is no attachment and detachment of STPs. This is no 
problem after all, because uia and uJd scale like with 
system size L and we are only looking at a local quan- 
tity (the shock position), so they can be neglected for 
large systems (which we are interested in). Besides, for 
real (finite) systems they are negligibly small compared 
to the other rates ijJh,oJs,ojf and ujb- On the other hand, 
if the STPs were allowed to detach, the undesirable pos- 
sibility of losing all the STPs through detachments could 
not be ruled out. Moreover, allowing STPs to attach and 
detach like the real particles would involve further sub- 
tleties of normalization during computation of averaged 
quantities. 

A STP, which is not located at the shock, has a ten- 
dency to move to the shock position. Moreover, if a STP 
is already located at the shock, it follows the shock as 
the shock moves. For the purpose of illustration, consider 
first an idealized shock of the form ...OOOOXXXXXX... In- 
serting a STP in cither the low density region or the high 
density region it is obvious from the rules given in (j54p 
that it will, on the average, move in the direction of the 
shock. However, in our model, the observed shocks are 
not ideal. Instead a few particles (holes) will appear in 
the low (high) density region. As a first approximation. 
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one can assume that these particles (holes) are isolated, 
e.g. configurations hke ...OOXOOOXXXOXX... Again, by 
careful use of the rules ([54| . one can show that the pre- 
ferred motion of the STP is towards the location of the 
shock also in such realistic situations. This argument can 
be refined even further. In appendix B we present an 
analytical argument in mean-field approximation which 
supports the heuristic arguments used in the illustrative 
examples in this paragraph. 

In addition to the rules listed above we define the fol- 
lowing fusion rules: 



2X 
2X 
X2 



01 with LOf 

02 with LOi, 
20 with tJb 



(55) 



The fusion rules ensure that if a shock exists, there will 
be a single STP in the system after sufficiently long time. 
This rule is extremely convenient because the lone STP 
will uniquely define the position of the sharp shock rather 
than a wide region of contiguous STPs separating the 
high-density and low-density regions. 

For the practical implementation of the STPs on the 
computer, one has to select the initial positions of the 
STPs. We chose to put one STP at each end of the sys- 
tem at the beginning of the simulation. If a shock can 
exist in the system, the STPs move to the shock position, 
fuse and, finally, indicate the shock position. We deter- 
mined the shock position in the stationary state by av- 
eraging over the fiuctuating positions of the lone STP in 
the steady state. In contrast, survival of two STPs in the 
steady state of the system indicates absence of any shock; 
instead, these two STPs indicate the formation of bound- 
ary layers. Although the latter phenomenon could be in- 
teresting, we shall not discuss it here. We have compared 
the shock position obtained following the STP approach 
with that inferred from the density profiles measured by 
computer simulations of our model. These comparisons 
established that the rules ([M)) and ([SS]). indeed, yield the 
correct results. 

We determined numerically the mean position of shock 
in a system with L = 600 sites as a function of Ua and 
Wfi which is shown as a 3D-plot in Fig. 1101 

In Fig. [Tl] we have plotted the shock positions as a 
function of the parameter c for different choices of lOs and 
LOf. In this figure, one observes two plateaus connected 
by a decaying domain (leftshift of the shock position). 
It seems that upper plateau approximately ends for c = 
l — Lds—LOf. Further detailed investigations will be needed 
to decide whether the sharp change in the position of the 
shock at this value of c indicates merely a crossover or a 
signature of a genuine phase transition. 



VII. ANALYTICAL PHASE DIAGRAM 
WITHOUT LANGMUIR DYNAMICS 

In this section we derive the phase diagram in the plane 
spanned by the boundary rates a and (3 for the special 





r- 600 

500 
400 
300 
200 
100 




i^a [1/s] 



FIG. 10: Shock position as function of Ua and Uh obtained 
from STP simulations. 



case of our model where attachments and detachment of 
the motors do not take place. In other words, we derive 
the phase diagram of our model in the a — /3-plane in 
the absence of Langmuir kinetics. We use the domain 
wall theory proposed in [s^l to derive this phase diagram 
from the flow-density relation (j26p of the corresponding 
periodic system. From this study, one can calculate the 
collective velocity and the shock velocity which determine 
the dynamics of the density profiles of the open system. 
Note that, because of the translational invariance of the 
periodic system, S and W show constant density profiles. 
The collective velocity Vc of this system is given by 



dJ{p) 

dp 



JfUJh 



CLOfP^ — 2[ljJh + UJs +^f)P + ^h+^s+^f 



Thus Vc — Q gives the critical density 



Pc = k - \/k{k - 1) 



where 



k = 



UJh +UJs+i^f 



(56) 



(57) 



(58) 



for the case c 7^ 0, and pc = ^/2 for the case c = 0. Note 
that k is always larger than 1. Next we calculate the 
shock velocity 



S = 



J{pl) - J{pr) 



(59) 



PL - PR 

where we take pL — a and pn = \ — (3. Then we have 

S^{ujh+ LOs + iOf){l3 -a) + cwfail - (3). (60) 
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FIG. 12: Phase diagram of this system without Langmuir 
dynamics. 



VIII. EXPERIMENTAL INVESTIGATION 
WITH KIFIA 




FIG. 11: Variation of the shock positions with the interaction 
parameter c for system size L — 3000, uid = 0.1 and: top: 
uOs = 55, ijJf = 145, LOa = 0.007, loh = 180; middle: ujs — 
100, ujf = 200, LOa = 0.1, LOh = 300; bottom: los = 100, Uf = 

300, LOa = 0.01, LOh = 130 



From S — we obtain the first order phase transition 
curve 



-k 



(61) 



that starts at a = and ends at a = Pc (Fig. HH) . This 
curve separates the low and high density phase. 



In the experiments performed by Okada [18| . micro- 
tubules labeled with a green fluorescent dye were immo- 
bilized on the top surface of the cell. The single-headed 
kinesins labeled with a red fluorescent dye were then in- 
troduced into the cell together with with ATP. The move- 
ment of the motor proteins was observed using imaging 
techniques of optical microscopy described in ^191] . A 
"comet-like" structure, as shown in Fig. [T21 was formed 
by the kinesins (red) on the microtubule (green). The 
first two images from the top, which correspond to low 
and moderate densities, respectively, were taken under 
essentially same conditions, but the lowermost image in 
the figure was taken with smaller intensifier gain, because 
it is too bright for the intensifier. 

No special filtering was applied to the original image. 
Each red fluorescent spot in Fig . [T51 normally corresponds 
to a single fluorescently-labeled kinesin molecule, if the 
density is not too high (top panel of Fig. [12] is a typi- 
cal example of such cases). Due to the optical resolu- 
tion limit (about 500 nm), more than one kinesin can 
together form a single brighter spot when the motors are 
too close to be resolved (as happens, for example, in the 
middle panel of Fig. [T3)) . Nevertheless, even in such sit- 
uations, the number of fluorochromes in each spot can 
be estimated from its intensity. At much higher densities 
(for example, that corresponding to the bottom panel of 
Fig. [13]), fluorescent signals are no longer separable as 
spots. Even in such cases, the density of fluorochromes 
can be estimated from their intensity profile. However, 
Okada measured the intensity profile just to confirm that 
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each spot corresponds to a single kinesin molecule in the 
lowest density experiment. In other words, at low densi- 
ties, the density of the fluorescent spots gives a good esti- 
mate for the density of kinesins. But, at higher densities, 
the spot density gives an underestimate of the kinesin 
density due to the overlap of fluorescent spots (which are 
not visually separable because of the limited resolution) . 

It is true that, under normal physiological conditions, 
the global density of motors in a cell never oversaturates 
the microtubule surface as happened in Okada's exper- 
iment described above. However, so far as the in-vivo 
situations are concerned, the motors and microtubules 
are heterogeneously distributed in cells. Thus, the local 
density of motors and microtubule surfaces might be a 
direct determinant of the formation of motor traffic jam 
within cells during in-vivo experiments. Moreover, in 
pathological situations, traffic jam on microtubule-based 
transport systems, such as axonal transport, is not rare. 
In fact, such traffic jams have been implicated in many 
neurodegenerative diseases [H, 113, HI]. Many putative 
factors may contribute to the "jammorigenesis" ; these 
include the population of the active motor proteins, the 
presence of the inactive motor proteins, the number of 
"obstacles" on the microtubule surface such as micro- 
tubule associated proteins, and so on. Obviously, these 
factors should be, ultimately, incorporated into a more 
"realistic" extended version of our model in order to ex- 
plicitly account for the observed "jammorigenesis" . The 
current version of our model is just the minimal one. 

These experimental results have three important impli- 
cations. First, traffic jam can actually take place in living 
cells at least in some experimental conditions. Second, 
the local concentration, rather than the global concentra- 
tion, of the motors determines whether or not jam will 
form in a living cell. Even in the overexpressing cells, 
the overall concentration of motors is much lower than 
that of tubulin. But still "comet" is formed. Third, neg- 
ative regulation systems, which are not included in the 
current version of our model, prevent jam formation in 
physiological situations. 



IX. DISCUSSION 

In this paper we have proposed a biologically moti- 
vated extension of our recent quantitative model de- 
scribing traffic-like collective movement of single-headed 
kinesin motors KIFIA. The dynamics of the system has 
been formulated in terms of a stochastic process where 
position of a motor is repesented by a discrete variable 
and time is continuous. The model explicitly captures 
the most essential features of the biochemical cycle of 
each motor by assigning two discrete internal ("chem- 
ical" or "conformational") states to each motor. The 
model not only takes into account the exclusion inter- 
actions, as in the previous models, but also includes a 
possible interaction of motors that controls ADP release 
rates by introducing a free parameter c. To our knowl- 




FIG. 13; (Color online) A "comet-like" structure formed by 
kinesins (red) on the microtubule (green). The high-density 
and low-density regions are clearly separated in this image. 
The white bar length is 2 /im. 

edge [27j , it is not possible even to establish the existence 
of this mechano-chemical interaction with the experimen- 
tal data currently available in the literature. However, we 
hope that our results reported here will help in develop- 
ing experimental methods which will not only test the 
existence of this interaction but also its strength if it ex- 
ists. For example, we have predicted the dependence of 
the shock position on c (and, therefore, that of c on the 
shock position). Thus, at least in prinicple, one could 
determine c by comparing the experimentally measured 
shock position with this relation. The c-dependence of 
some of the other quantities reported here may provide 
alternative, and possibly, more direct way of estimating 
the strength of this mechano-chemical interaction. 

We have compared and contrasted our model and the 
results with earlier generic models of single motors as well 
as those of motor traffic. Our analytical treatment of the 
dynamical equations in the continuum limit (i.e., a limit 
in which the spatial position of each motor is denoted by 
a continuous variable) has also established the occurrence 
of a non-propagating shock in this model. We have also 
calculated the position of this shock numerically using 
the method of second class particles. 

Mean field treatment of the rate equations for c = 
showed that this special case of our model is equivalent to 
the simpler PFF model which also predicts two-phase co- 
existence (where the two phases are separated by a non- 
propagating shock). One can argue analytically [H, H^l, 
as we also observed in simulations, that the general fea- 
tures of the a — /3-phase diagrams of our model is the 
same as those for the PFF model. Thus, the PFF model, 
in spite of its simplicity, captures the essential generic 
features of intracellular transport. But, it is not possi- 
ble to make direct quantitative comparison between the 
predictions of the PFF model and experimental data as 
the parameters of the PFF model are not accessible to 
direct biochemical experiments. In contrast, our model 
captures the essential features of the internal biochemical 
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transitions of each single-headed kinesin and we could es- 
tablish a one to one correspondence between our model 
parameters and measurable quantities. The concentra- 
tions of the kinesin motors and ATP are two such impor- 
tant parameters both of which which are variable in-vivo 
and can be controlled in m-i)itro-experiments. We have 
reported the phase diagram of our model in the plane 
spanned by these experimentally accessible parameters. 

Finally, we have summarized evidences for the forma- 
tion of molecular motor jam from Okada's in- vitro exper- 
iments [18] and discussed their relevance in intra-cellular 
transport under physiological conditions. 
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Appendix A: A CENERALIZED MODEL OF NON- 
INTERACTING MOTORS 

In order to make a comparison between the non- 
interacting limit of our model and the earlier models 
of non-interacting molecular motors, we consider here a 
slightly more general model which allows "reverse" tran- 
sitions for each of the "forward" transitions. Then we 
show that the non-interacting limit of our model is a 
special case of the general model while some other spe- 
cial cases correspond to earlier models of non-interacting 
motors. 

Consider the multi-step chemical kinetic scheme shown 
in the Fig. [TH Note that this generalized scheme [11] al- 
lows a transition from the strongly bound state at i -I- 1 
to the weakly bound state at i with the rate constant ujj 
which is not allowed in our model shown in fig (21 In fact, 
in this generalized scheme, corresponding to every for- 
ward step (those corresponding to ujj', ui'^ and uib) there 
is a backward step (corresponding to tuj, and uJb, re- 
spectively). This generalization is in the spirit of Fisher- 
Kolomeisky-type multi-step chemical kinetic models of 
molecular motors [H, [s^, [sij where each of the reactions 
are allowed to be reversible, albeit with different rate 
constants, in general. 

In the mean-field limit the the master equations gov- 
erning the dynamics of this general model in the bulk are 



FIG. 14: Schematic description of a general three-state model 
of a single molecular motor. 



given by 



dt 



dW,, 



+ LOf^Wi - uj^St - ujjSi - LOdSi, (62) 



+ LUbiW,-i+W,+i)-2LUbW,. (63) 

Imposing periodic boundary conditions, the steady state 
solutions for S and W can be written as 



S 



W 



U)a{i^h + + + ) + Wrf(tJ^ + Ujj) 



(64) 



(65) 



Hence, 



S + W= —— y ,f f- XT .(66) 

CJq(cJ^ + UjJ + UJ]^ + UJj ) + UJd[tJ^ + LjJ) 

The corresponding steady-state flux 



is given by 

J = 



{K + l){uJi,+uj+) + [w++^J) 



(67) 



(68) 



The relation between this generalized model of 
non-interacting motors and the non-interacting limit of 
our model is quite straightforward. In the special case 
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— oj/i and 



= 0, using the identification wj" = ujf, uj^^ 
uj^ = uJs, the equations (IM)) , and reduce to the 
equations pT)) . p8| and (flQl) . respectively. 

Appendix B: A MF ARGUMENT FOR MOVEMENT 
OF STP AND SHOCK POSITION 



where we have used the fact that 

i 

= iuf VF(i)(^I¥,+i-^W,_i) =0. (74) 



In this appendix we argue that a STP will move to the 
location of the shock, if a shock exists in the system. Our 
arguments are based on an analysis in the mean-field ap- 
proximation. The master equations for the probabilities 
of the STP, which correspond to the equations ([9]) for the 
real particles, are given by 



dt 



- u^sf + u^M^^ + ^/(l - c)M^f VS, (69) 



— W 



(2) 



-c.Kl-c)VFfV!i\, 



.(2) 



(70) 



(2) (2) 

where SI and represent the probabilities of finding 
the STP in the weakly (W^^)) and strongly (S'^^)) bound 

states, respectively, at the site i] note that sf''' and wl^'^ 
are the corresponding probabilities for the real particles. 

Obviously, = S-f '^^ + W'l^''^\ 

Adding the two equations (|69p and ((70)l we obtain 



dt 



(2) , Ti/(2) /I Ti^(l) (2)\ 

- (^;M^f)(l-pW)-c./MA«pl+\). (71) 

Comparing this with the equation of continuity dpi j dt — 
Ji_i — Ji,^e identify the current of STP on site i 

to be 



J, 



STP 



iOfWPil-f^^-cofWl'^p'^,. (72) 

Consider a situation where we have one STP in a con- 
tinuous region of particles (with no shock inside), so we 

(1) ^ St\ =: SW and 



can put p. 



(1) 



(1) 



wl^\ =: W^^'). We assume that, after sufii- 



Ji) Jl). Q 



ciently long time, the internal states of the STP relax to 
a stationary state so that the probabilities of finding the 
STP in the strongly-bound and weakly-bound states are 
independent of time. However, the mean position of the 
STP might still change with time. 

Then, S := Si is the probability of finding the 
STP in a strongly bound state, while the corresponding 
probability of finding the STP in the weakly bound state 
is W := Wj^ where the summations are over an in- 
terval of length I that contains no shocks and one single 
STP. Obviously, S + W = 1. Using ^ we have 



w/(l-c)p(i)W, (73) 



Solving equation ([75)1 for W, we obtain 



W = 



Cl)/i +UJs + {I - Cp'^'^'>)uJf 



(75) 



The results derived above are valid for any density dis- 
tribution of STPs as long as there is a shock-free neigh- 
bourhood of the STP and the particles are in a steady 
state. Now consider a specific configuration where a STP 
is given to be located at the site i while its internal state 
remains unspecified. In this case, 



(2) X 
Pk = Oik- 



(76) 



(2) 

Then we have W = for any summation interval that 

('21 

includes the site i. Of course, Wl^ ' = for fc 7^ i for this 
distribution of pk- Therefore, using (j75p we obtain 



W, 



(2) 



The analogous solution for VF'^^ obtained from (|27p is 



(77) 



UJhP 



(1) 



UJh + UJs + {I - Cp^^'>)LOf 



(78) 



For the density distribution considered here, we can 
take the current as an effective hopping rate of the STP 

f 21 

to the right, i.e., (j[ = Ji{{Pk ^ik})- Similarly, we 
have effective hopping rate of the STP to the left q\ — 
-Ji-i{{Pk'' = Sik})- Inserting ^ and ^ into ^ 
for = (5jfc, we obtain 



JfUJh 



UJh+ UJs + {'^ ~ Cp^'^'))ujf 



(l-2pW). (79) 



Note that the fraction in ([79|l is always positive. There- 
fore, if a STP is in a low density region with p^^^ < ^, 
we have q[ — > and the STP tends to hop to the 
right. But, if the STP is in a high density region with 
/d^^-* > i, we have — ql < and its prcfcrcd direction 
of hopping is left. Thus, in the continuum limit, if there 
is one shock separating a low density region at the left 
and a high density region at the right, any single STP 
will be driven to this domain wall. For sufficiently long 
time the average position of the STP will be equal to the 
shock position. 
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